classdef EddyDatasetClass < handle
    % Eddy 实例类 储存Eddy数据
    % 使用导入 EddyFile 文件初始化 

    properties
        EddyFile
        Dict
        traj
        % 特别注意这个extent是当前操作的extent，可以通过改变extent再进行func变换研究区域
        extent = [-90 90;-180 180];
        year = [1993  2020];
    end
    
    methods
        % 通过导入涡旋文件初始化构造函数
        function obj = EddyDatasetClass(EddyFile, varargin)
            % 构造函数 初始化类 varargin 包括：extent year
            obj.EddyFile = EddyFile;               
            idx_extent = find(strcmp(varargin,'extent'));
            if ~isempty(idx_extent)
                obj.extent = varargin{idx_extent + 1};
            end

            idx_year = find(strcmp(varargin,'year'));
            if ~isempty(idx_year)
                obj.year = varargin{idx_year + 1};
            end
            % obj.eddyFile2Dict();
            % obj.Dict2traj();
        end

        % 设置extent
        function extentSet(obj, extent)
            obj.extent = extent;
        end

        % 设置year
        function yearSet(obj, year)
            obj.year = year;
        end

        % % 从file中读取eddy转换呈dict
        function eddyFile2Dict(obj)
            obj.Dict = EddyDatasetClass.setDict();
            latitude = ncread(obj.EddyFile,'latitude');
            longitude = ncread(obj.EddyFile,'longitude');
            track = ncread(obj.EddyFile,'track');
            observation_number = ncread(obj.EddyFile,'observation_number');
            Time = ncread(obj.EddyFile,'time');


            selectIndex = Time >= days(datetime(obj.year(1),1,1) - datetime(1950,1,1)) & Time < days(datetime(obj.year(2) + 1,1,1) - datetime(1950,1,1));
            selectTrack = unique(track(selectIndex));
            index = find(ismember(track, selectTrack) == 1);
            fields = fieldnames(obj.Dict);
            
            for i = 1:length(fields)
                field = fields(i);
                key = field{1};
                data = ncread(obj.EddyFile,key);
                [m,n] = size(data);
                if n == 1
                    if strcmp(field,'time')
                        obj.Dict.(key) = datetime(1950,1,1) + days(data(index));
                    else
                        obj.Dict.(key) = data(index);
                    end   
                else
                    obj.Dict.(key) = data(1:m,index);
                end
            end

            obj.Dict2traj();
            idx = [];
            for i = 1:length(obj.traj)
                i
                traj_item = obj.traj(i);
                longitude = traj_item.longitude;
                latitude = traj_item.latitude;
                %%%加一层振幅和生命周期过滤
                lifetime = length(longitude);
                %mean_amplitude = mean(traj_item.amplitude);
                if lifetime <28 
                    continue;
                end
                %%% 适应不同的extent 额外添加对[-90 90;0 360]这种格式的extent的支持
                [~,n] = size(obj.extent);
                if n == 2
                    xboundary = obj.extent(2,:);
                    yboundary = obj.extent(1,:);
                    xbox = xboundary([1 1 2 2 1]);
                    ybox = yboundary([1 2 2 1 1]);
                    in = inpolygon(longitude,latitude,xbox, ybox);
                else
                    in = inpolygon(longitude,latitude,obj.extent(2,:), obj.extent(1,:));
                end
                %%%
                if all(in(:)==1)
                    idx = [idx i];
                end
            end
            obj.traj = obj.traj(idx);
            obj.Traj2dict();

        end
        

        %将Dict 转换成 traj（涡旋轨迹）支持extent
        function Dict2traj(obj)
            %读取Dict
            LatitudeArray = obj.Dict.latitude;
            LongitudeArray = obj.Dict.longitude;
            TrackArray = obj.Dict.track;
            unique_track = unique(TrackArray);
            for i = 1:length(unique_track)
                trajectory_Index = find(obj.Dict.track==unique_track(i));
                [traj_item] = eddyDictByIndex(obj.Dict, trajectory_Index);  
                traj_item.lifetime = length(traj_item.time);
                trajectory(i) = traj_item;
            end
            obj.traj = trajectory;
            % 函数根据index 切割Dict
            function [newDict] = eddyDictByIndex(oldDict, index)
                fields = fieldnames(oldDict);
                for fieldI = 1:length(fields)
                    field = fields(fieldI);
                    key = field{1};
                    data = oldDict.(key);
                    [m,n] = size(data);
                    if n == 1
                        newDict.(key) = data(index);
                    else 
                        newDict.(key) = data(1:m,index);
                    end
                end
            end  
        end

        function Traj2dict(obj)
            dict = obj.traj(1);
            for i = 2:length(obj.traj) 
                i
                traj_item = obj.traj(i);
                dict = DictAdd(dict, traj_item);

            end
            obj.Dict = dict;
            function [newDict] = DictAdd(Dict, traj_item)
                fields = fieldnames(Dict);
                for fieldI = 1:length(fields)
                    field = fields(fieldI);
                    if strcmp(field,'lifetime')
                        continue;
                    end
                    key = field{1};
                    adddata = traj_item.(key);
                    data = Dict.(key);
                    [m,n] = size(data);
                    if n == 1
                        data = [data; adddata];
                        newDict.(key) = data;
                    else 
                        data = [data, adddata];
                        newDict.(key) = data;
                    end
                end
            end
        end
    end




    methods(Static)
        %创建一个涡旋字典
        function dict = setDict()
            % 初始化构造字典
            dict = struct('amplitude',0,...
            'cost_association',0,...
            'effective_area',0,...
            'effective_contour_height',0,...
            'effective_contour_latitude',0,...
            'effective_contour_longitude',0,...
            'effective_contour_shape_error',0,...
            'effective_radius',0,...
            'inner_contour_height',0,...
            'latitude',0,...
            'latitude_max',0,...
            'longitude',0,...
            'longitude_max',0,...
            'num_contours',0,...
            'num_point_e',0,...
            'num_point_s',0,...
            'observation_flag',0,...
            'observation_number',0,...
            'speed_area',0,...
            'speed_average',0,...
            'speed_contour_height',0,...
            'speed_contour_latitude',0,...
            'speed_contour_longitude',0,...
            'speed_contour_shape_error',0,...
            'speed_radius',0,...
            'time',0,...
            'track',0,...
            'uavg_profile',0);
        end

        % Dict 转成 有独立涡旋形成的数组
        function eddyArray =  eddyArrayAddByDict(eddyArray, Dict, type)
            Array_Len = length(eddyArray);
            for i = 1:length(Dict.time)
                eddyArray(Array_Len + i).latitude = Dict.latitude(i);
                eddyArray(Array_Len + i).longitude = Dict.longitude(i);
                eddyArray(Array_Len + i).contour_latitude = Dict.effective_contour_latitude(:,i);
                eddyArray(Array_Len + i).contour_longitude = Dict.effective_contour_longitude(:,i);
                eddyArray(Array_Len + i).radius = Dict.effective_radius(i);
                eddyArray(Array_Len + i).time = Dict.time(i);
                eddyArray(Array_Len + i).type = type;
                eddyArray(Array_Len + i).amplitude = Dict.amplitude(i);
                % eddyArray(Array_Len + i).speed = trajectoryItem.speed_average(i);
            end
        end
        
        % 将每涡旋轨迹转换成由每个独立涡旋形成的数组
        function eddyArray = eddyArrayAddByTraj(eddyArray, trajectory, type)
            for i = 1:length(trajectory)
                n = length(eddyArray);
                trajectoryItem = trajectory(i);
                [len , ~]= size(trajectoryItem.time);
                % Begin = floor(len / 5);
                % End = floor(4 * len / 5);
                count = 0;
                % for j = Begin:End
                for j = 1:len
                    count = count + 1;
                    eddyArray(n + count).latitude = trajectoryItem.latitude(j);
                    eddyArray(n + count).longitude = trajectoryItem.longitude(j);
                    eddyArray(n + count).contour_latitude = trajectoryItem.effective_contour_latitude(:,j);
                    eddyArray(n + count).contour_longitude = trajectoryItem.effective_contour_longitude(:,j);
                    eddyArray(n + count).radius = trajectoryItem.effective_radius(j);
                    eddyArray(n + count).time = trajectoryItem.time(j);
                    eddyArray(n + count).type = type;
                    eddyArray(n + count).amplitude = trajectoryItem.amplitude(j);
                    eddyArray(n + count).speed = trajectoryItem.speed_average(j);
                end
            end
        end

                %%% 将轨迹按照季节拆分
        % 拆分矩阵格式 [12 1 2;3 4 5;6 7 8;9 10 11]
        function [spilted_traj] = trajSplitBySeason(traj, Traj_Cluster_Matrix)
            [m,n] =size(Traj_Cluster_Matrix);
            % 分割结果矩阵进行置空
            for i = 1:m
                spilted_traj(i).traj = [];
            end
            
            for i = 1:length(traj)
                traj_Item = traj(i);
                % 计算轨迹的长度
                taij_len = length(traj_Item.time);

                % 定义两个指针 一个起始指针 一个结束指针
                begin_cursor = 1;
                end_cursor = 1;

                % 设置循环 只要结束指针不到末尾 就不断切割traj
                while end_cursor < taij_len
                    end_cursor = end_cursor + 1;
                    begin_Cluster_ID = judgDateInCluster(traj_Item.time(begin_cursor), Traj_Cluster_Matrix);
                    end_Cluster_ID = judgDateInCluster(traj_Item.time(end_cursor), Traj_Cluster_Matrix);
                    
                    if begin_Cluster_ID ~= end_Cluster_ID
                        [spiled_traj_item] = eddyDictByIndex(traj_Item, begin_cursor:end_cursor - 1);  
                        spiled_traj_item.lifetime = length(spiled_traj_item.time);
                        spilted_traj(begin_Cluster_ID).traj = [spilted_traj(begin_Cluster_ID).traj spiled_traj_item];
                        begin_cursor = end_cursor;
                    end     
                    if end_cursor == taij_len
                        [spiled_traj_item] = eddyDictByIndex(traj_Item, begin_cursor:end_cursor);  
                        spilted_traj(begin_Cluster_ID).traj = [spilted_traj(begin_Cluster_ID).traj spiled_traj_item];
                        begin_cursor = end_cursor;
                    end
                end
                    
            end
            % 判断涡旋日期在哪个Cluster聚类里，输出聚类的序号
            function [Cluster_ID] = judgDateInCluster(date, Traj_Cluster_Matrix)
                [m, ~] = size(Traj_Cluster_Matrix);
                for p = 1:m
                    month = date.Month;
                    if ismember(month, Traj_Cluster_Matrix(p,:))
                        Cluster_ID = p;
                        return
                    end
                end
            end
            % 函数根据index 切割Dict
            function [newDict] = eddyDictByIndex(oldDict, index)
                fields = fieldnames(oldDict);
                for fieldI = 1:length(fields)
                    field = fields(fieldI);
                    key = field{1};
                    data = oldDict.(key);
                    [m,n] = size(data);
                    if n == 1 & m == 1
                        newDict.(key) = data;
                        
                    elseif n == 1 & m ~= 1
                        newDict.(key) = data(index);
                    else
                        newDict.(key) = data(1:m,index);
                    end
                end
            end   
        end
    end
end

